home *** CD-ROM | disk | FTP | other *** search
/ QRZ! Ham Radio 8 / QRZ Ham Radio Callsign Database - Volume 8.iso / pc / files / ant_nec / nec81tar.z / nec81tar / fbngf.f < prev    next >
Text File  |  1991-05-13  |  22KB  |  775 lines
  1. C $TITLE: 'FBNGF'
  2. C $NOFLOATCALLS
  3.       SUBROUTINE FBNGF (NEQ,NEQ2,IRESRV,IB11,IC11,ID11,IX11,IW)
  4. C     FBNGF SETS THE BLOCKING PARAMETERS FOR THE B, C, AND D ARRAYS
  5. C     FOR OUT-OF-CORE STORAGE.
  6.       INTEGER*4 NEQ,NEQ2,IB11,IC11,ID11,IX11,NBLN,NDLN,NBCD,IR,IRESX
  7.       INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  8.       COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
  9.      1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  10. C**
  11. C     D      WRITE(*,*) '  FBNGF: START'
  12. C**
  13.       IRESX=IRESRV-IMAT
  14.       NBLN=NEQ*NEQ2
  15.       NDLN=NEQ2*NEQ2
  16.       NBCD=2*NBLN+NDLN
  17.       IF (NBCD.GT.IRESX) GO TO 1
  18.       ICASX=1
  19.       IB11=IMAT+1
  20.       GO TO 2
  21. 1      CONTINUE
  22. C**
  23. C     D      WRITE(*,*) '  FBNGF: OPEN 11,12,13,14,15'
  24. C**
  25.       OPEN (11,FORM='UNFORMATTED')
  26.       OPEN (12,FORM='UNFORMATTED')
  27.       OPEN (13,FORM='UNFORMATTED')
  28.       OPEN (14,FORM='UNFORMATTED')
  29.       OPEN (15,FORM='UNFORMATTED')
  30.       IF (ICASE.LT.3) GO TO 3
  31.       IF (NBCD.GT.IRESRV.OR.NBLN.GT.IRESX) GO TO 3
  32.       ICASX=2
  33.       IB11=1
  34. 2     NBBX=1
  35.       NPBX=NEQ
  36.       NLBX=NEQ
  37.       NBBL=1
  38.       NPBL=NEQ2
  39.       NLBL=NEQ2
  40.       GO TO 5
  41. 3     IR=IRESRV
  42.       IF (ICASE.LT.3) IR=IRESX
  43.       ICASX=3
  44. C**
  45. C     D      WRITE(*,*) '  FBNGF: OPEN 16'
  46. C**
  47.       OPEN (16,FORM='UNFORMATTED')
  48.       IF (NDLN.GT.IR) ICASX=4
  49.       NBCD=2*NEQ+NEQ2
  50.       NPBL=IR/NBCD
  51.       NLBL=IR/(2*NEQ2)
  52.       IF (NLBL.LT.NPBL) NPBL=NLBL
  53.       IF (ICASE.LT.3) GO TO 4
  54.       NLBL=IRESX/NEQ
  55.       IF (NLBL.LT.NPBL) NPBL=NLBL
  56. 4     IF (NPBL.LT.1) GO TO 6
  57.       NBBL=(NEQ2-1)/NPBL
  58.       NLBL=NEQ2-NBBL*NPBL
  59.       NBBL=NBBL+1
  60.       NBLN=NEQ*NPBL
  61.       IR=IR-NBLN
  62.       NPBX=IR/NEQ2
  63.       IF (NPBX.GT.NEQ) NPBX=NEQ
  64.       NBBX=(NEQ-1)/NPBX
  65.       NLBX=NEQ-NBBX*NPBX
  66.       NBBX=NBBX+1
  67.       IB11=1
  68.       IF (ICASE.LT.3) IB11=IMAT+1
  69. 5     IC11=IB11+NBLN
  70.       ID11=IC11+NBLN
  71.       IX11=IMAT+1
  72.       WRITE(IW,11)  NEQ2
  73.       IF (ICASX.EQ.1) RETURN
  74.       WRITE(IW,8)  ICASX
  75.       WRITE(IW,9)  NBBX,NPBX,NLBX
  76.       WRITE(IW,10)  NBBL,NPBL,NLBL
  77. C**
  78. C     D      WRITE(*,*) '  FBNGF: RETURN AT END'
  79. C**
  80.       RETURN
  81. 6     WRITE(IW,7)  IRESRV,IMAT,NEQ,NEQ2
  82.       WRITE(*,7)  IRESRV,IMAT,NEQ,NEQ2
  83.       STOP
  84. C
  85. 7     FORMAT(' ERROR - INSUFFICIENT STORAGE FOR INTERACTION MATRICIES',
  86.      1'  IRESRV,IMAT,NEQ,NEQ2 =',4I5)
  87. 8     FORMAT(' FILE STORAGE FOR NEW MATRIX SECTIONS -  ICASX =',I2)
  88. 9     FORMAT(' B FILLED BY ROWS-',15X,'NO. BLOCKS =',I3,3X,
  89.      1'ROWS PER BLOCK =',I3,3X,'ROWS IN LAST BLOCK =',I3)
  90. 10    FORMAT(' B BY COLUMNS, C AND D BY ROWS -',2X,'NO. BLOCKS =',
  91.      1 I3,4X,'R/C PER BLOCK =',I3,4X,'R/C IN LAST BLOCK =',I3)
  92. 11    FORMAT (//,' N.G.F. - NUMBER OF NEW UNKNOWNS IS',I4)
  93.       END
  94. C
  95. C
  96. C
  97.       SUBROUTINE FBLOCK (NROW,NCOL,IMAX,IRNGF,IPSYM,IW)
  98. C     FBLOCK SETS PARAMETERS FOR OUT-OF-CORE SOLUTION FOR THE PRIMARY
  99. C     MATRIX (A)
  100.       INTEGER*4 IMX1,NROW,NCOL,IPSYM,KK,I,J,K,KA,NOP
  101.       INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  102.       REAL*8 ARG,PHAZ
  103.       COMPLEX*16 SSX,DETER
  104.       COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,
  105.      1 IMAT,ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  106.       COMMON /SMAT/ SSX(16,16)
  107. C     D      WRITE(*,*) '  FBLOCK: START'
  108.       IMX1=IMAX-IRNGF
  109.       IF (NROW*NCOL.GT.IMX1) GO TO 2
  110.       NBLOKS=1
  111.       NPBLK=NROW
  112.       NLAST=NROW
  113.       IMAT=NROW*NCOL
  114.       IF (NROW.NE.NCOL) GO TO 1
  115.       ICASE=1
  116. C     D      WRITE(*,*) '  FBLOCK: RETURN BEFORE 1'
  117.       RETURN
  118. 1     ICASE=2
  119.       GO TO 5
  120. 2      CONTINUE
  121. C**
  122. C     D      WRITE(*,*) '  FBLOCK: OPEN 11,12,13,14'
  123. C**
  124.       OPEN (11,FORM='UNFORMATTED')
  125.       OPEN (12,FORM='UNFORMATTED')
  126.       OPEN (13,FORM='UNFORMATTED')
  127.       OPEN (14,FORM='UNFORMATTED')
  128.       IF (NROW.NE.NCOL) GO TO 3
  129.       ICASE=3
  130.       NPBLK=IMAX/(2*NCOL)
  131.       NPSYM=IMX1/NCOL
  132.       IF (NPSYM.LT.NPBLK) NPBLK=NPSYM
  133.       IF (NPBLK.LT.1) GO TO 12
  134.       NBLOKS=(NROW-1)/NPBLK
  135.       NLAST=NROW-NBLOKS*NPBLK
  136.       NBLOKS=NBLOKS+1
  137.       NBLSYM=NBLOKS
  138.       NPSYM=NPBLK
  139.       NLSYM=NLAST
  140.       IMAT=NPBLK*NCOL
  141.       WRITE(IW,14)  NBLOKS,NPBLK,NLAST
  142.       GO TO 11
  143. 3     NPBLK=IMAX/NCOL
  144.       IF (NPBLK.LT.1) GO TO 12
  145.       IF (NPBLK.GT.NROW) NPBLK=NROW
  146.       NBLOKS=(NROW-1)/NPBLK
  147.       NLAST=NROW-NBLOKS*NPBLK
  148.       NBLOKS=NBLOKS+1
  149.       WRITE(IW,14)  NBLOKS,NPBLK,NLAST
  150.       IF (NROW*NROW.GT.IMX1) GO TO 4
  151.       ICASE=4
  152.       NBLSYM=1
  153.       NPSYM=NROW
  154.       NLSYM=NROW
  155.       IMAT=NROW*NROW
  156.       WRITE(IW,15)
  157.       GO TO 5
  158. 4     ICASE=5
  159.       NPSYM=IMAX/(2*NROW)
  160.       NBLSYM=IMX1/NROW
  161.       IF (NBLSYM.LT.NPSYM) NPSYM=NBLSYM
  162.       IF (NPSYM.LT.1) GO TO 12
  163.       NBLSYM=(NROW-1)/NPSYM
  164.       NLSYM=NROW-NBLSYM*NPSYM
  165.       NBLSYM=NBLSYM+1
  166.       WRITE(IW,16)  NBLSYM,NPSYM,NLSYM
  167.       IMAT=NPSYM*NROW
  168. 5      CONTINUE
  169.       NOP=NCOL/NROW
  170.       IF (NOP*NROW.NE.NCOL) GO TO 13
  171.       IF (IPSYM.GT.0) GO TO 7
  172. C
  173. C     SET UP SSX MATRIX FOR ROTATIONAL SYMMETRY.
  174. C
  175.       PHAZ=6.2831853072D0/NOP
  176.       DO 6 I=2,NOP
  177.       DO 6 J=I,NOP
  178.       ARG=PHAZ*(I-1)*(J-1)
  179.       SSX(I,J)=DCMPLX(DCOS(ARG),DSIN(ARG))
  180.       SSX(J,I)=SSX(I,J)
  181. 6      CONTINUE
  182.       GO TO 11
  183. C
  184. C     SET UP SSX MATRIX FOR PLANE SYMMETRY
  185. C
  186. 7     KK=1
  187.       SSX(1,1)=DCMPLX(1.,0.)
  188.       IF ((NOP.EQ.2).OR.(NOP.EQ.4).OR.(NOP.EQ.8)) GO TO 8
  189.       STOP
  190. 8     KA=NOP/2
  191.       IF (NOP.EQ.8) KA=3
  192.       DO 10 K=1,KA
  193.       DO 9 I=1,KK
  194.       DO 9 J=1,KK
  195.       DETER=SSX(I,J)
  196.       SSX(I,J+KK)=DETER
  197.       SSX(I+KK,J+KK)=-DETER
  198. 9     SSX(I+KK,J)=DETER
  199. 10    KK=KK*2
  200. 11      CONTINUE
  201. C**
  202. C     D      WRITE(*,*) '  FBLOCK: RETURN AFTER 11'
  203. C**
  204.       RETURN
  205. 12    WRITE(IW,17)  NROW,NCOL
  206.       WRITE(*,17)  NROW,NCOL
  207.       STOP
  208. 13    WRITE(IW,18)  NROW,NCOL
  209.       WRITE(*,18)  NROW,NCOL
  210.       STOP
  211. C
  212. 14    FORMAT (//' MATRIX FILE STORAGE-NO. BLOCKS=',I5,
  213.      1' COLUMNS PER BLOCK=',I5,' COLUMNS IN LAST BLOCK=',I5)
  214. 15    FORMAT(' SUBMATRICIES FIT IN CORE')
  215. 16    FORMAT(' SUBMATRIX PARTITIONING-NO. BLOCKS=',I5,
  216.      1' COLUMNS PER BLOCK=',I5,' COLUMNS IN LAST BLOCK=',I5)
  217. 17    FORMAT(' ERROR - INSUFFICIENT STORAGE FOR MATRIX',2I5)
  218. 18    FORMAT(' SYMMETRY ERROR - NROW,NCOL=',2I5)
  219.       END
  220. C
  221. C
  222. C
  223.       SUBROUTINE LOAD(ZARRAY,ZLR,ZLI,ZLC,SI,BI,LD,ITAG,LDTYP,
  224.      1 LDTAG,LDTAGF,LDTAGT,IW)
  225. C
  226. C     LOAD CALCULATES THE IMPEDANCE OF SPECIFIED SEGMENTS FOR VARIOUS
  227. C     TYPES OF LOADING
  228. C
  229. C***
  230. C***  FIX TO ADD PERMEABILITY TO WIRES - 12 DEC 88 - R W ADLER
  231. C***  ADD VALUE OF REL. PERM. AFTER CONDUCTIVITY ON LD CARD
  232. C***  2 LINES ADDED - 1 LINE CHANGED - ALSO FIXES TO FUNCTION ZINT
  233. C***
  234.       CHARACTER ADUM*20
  235.       REAL*8 TPCJX,ZREAL,ZIMAG
  236.       COMPLEX*16 TPCJ
  237.       COMPLEX*16 ZARRAY,ZT,ZINT
  238.       INTEGER*4 ITAG,N1,N2,N,NP,M1,M2,M,MP,IPSYM
  239.       COMMON/DATAD/N1,N2,N,NP,M1,M2,M,MP,IPSYM,WLAM
  240.       COMMON /ZLOAD/ NLOAD,NLODF
  241.       DIMENSION LDTYP(30),LDTAG(30),LDTAGF(30),LDTAGT(30),ZLR(30),
  242.      1 ZLI(30),ZLC(30),TPCJX(2),ZARRAY(LD),ITAG(LD),SI(LD),BI(LD)
  243.       EQUIVALENCE (TPCJ,TPCJX)
  244.       DATA TPCJX/0.,1.883698955D+9/
  245. C**
  246. C     D      WRITE(*,*) '  LOAD: START, NLOAD=',NLOAD,' NLODF=',NLODF
  247.       FZERO=0.
  248. C**
  249. C
  250. C     WRITE(IW,HEADING)
  251. C
  252.       WRITE(IW,25)
  253. C
  254. C     INITIALIZE D ARRAY, USED FOR TEMPORARY STORAGE OF LOADING
  255. C     INFORMATION.
  256. C
  257.       DO 1 I=N2,N
  258.       ZARRAY(I)=DCMPLX(0.,0.)
  259. 1      CONTINUE
  260.       IWARN=0
  261. C
  262. C     CYCLE OVER LOADING CARDS
  263. C
  264.       ISTEP=0
  265.  2    ISTEP=ISTEP+1
  266.       IF (ISTEP.LE.NLOAD) GO TO 5
  267.       IF (IWARN.EQ.1) WRITE(IW,26)
  268.       IF (N1+2*M1.GT.0) GO TO 4
  269.       NOP=N/NP
  270.       IF (NOP.EQ.1) GO TO 4
  271.       DO 3 I=1,NP
  272.       ZT=ZARRAY(I)
  273.       L1=I
  274.       DO 3 L2=2,NOP
  275.       L1=L1+NP
  276.       ZARRAY(L1)=ZT
  277. 3      CONTINUE
  278.  4    CONTINUE
  279. C**
  280. C     D      WRITE(*,*) '  LOAD: RETURN AFTER 4'
  281. C**
  282.       RETURN
  283.  5    IF (LDTYP(ISTEP).LE.5) GO TO 6
  284.       WRITE(IW,27)  LDTYP(ISTEP)
  285.       STOP
  286.  6    LDTAGS=LDTAG(ISTEP)
  287.       JUMP=LDTYP(ISTEP)+1
  288.       ICHK=0
  289. C
  290. C     SEARCH SEGMENTS FOR PROPER ITAGS
  291. C
  292.       L1=N2
  293.       L2=N
  294.       IF (LDTAGS.NE.0) GO TO 7
  295.       IF (LDTAGF(ISTEP).EQ.0.AND.LDTAGT(ISTEP).EQ.0) GO TO 7
  296.       L1=LDTAGF(ISTEP)
  297.       L2=LDTAGT(ISTEP)
  298.       IF (L1.GT.N1) GO TO 7
  299.       WRITE(IW,29)
  300.       STOP
  301. 7      CONTINUE
  302.       DO 17 I=L1,L2
  303.       IF (LDTAGS.EQ.0) GO TO 8
  304.       IF (LDTAGS.NE.ITAG(I)) GO TO 17
  305.       IF (LDTAGF(ISTEP).EQ.0) GO TO 8
  306.       ICHK=ICHK+1
  307.       IF (ICHK.GE.LDTAGF(ISTEP).AND.ICHK.LE.LDTAGT(ISTEP)) GO TO 9
  308.       GO TO 17
  309.  8    ICHK=1
  310. C
  311. C     CALCULATION OF LAMDA*IMPED. PER UNIT LENGTH, JUMP TO APPROPRIATE
  312. C     SECTION FOR LOADING TYPE
  313. C
  314. 9      CONTINUE
  315. C**
  316. C     D      WRITE(*,*) '  LOAD: AFTER  9, JUMP=',JUMP
  317. C**
  318.       GO TO (10,11,12,13,14,15), JUMP
  319. 10      CONTINUE
  320.       ZT=ZLR(ISTEP)/SI(I)+TPCJ*ZLI(ISTEP)/(SI(I)*WLAM)
  321.       IF(ABS(ZLC(ISTEP)).GT.1.E-20) ZT=ZT+WLAM/(TPCJ*SI(I)*ZLC(ISTEP))
  322.       GO TO 16
  323.  11   ZT=TPCJ*SI(I)*ZLC(ISTEP)/WLAM
  324.       IF(ABS(ZLI(ISTEP)).GT.1.E-20) ZT=ZT+SI(I)*WLAM/(TPCJ*ZLI(ISTEP))
  325.       IF (ABS(ZLR(ISTEP)).GT.1.E-20) ZT=ZT+SI(I)/ZLR(ISTEP)
  326.       ZT=1./ZT
  327.       GO TO 16
  328.  12   ZT=ZLR(ISTEP)*WLAM+TPCJ*ZLI(ISTEP)
  329.       IF (ABS(ZLC(ISTEP)).GT.1.E-20) ZT=ZT+1./(TPCJ*SI(I)*SI(I)*
  330.      1 ZLC(ISTEP))
  331.       GO TO 16
  332. 13      ZT=TPCJ*SI(I)*SI(I)*ZLC(ISTEP)
  333.       IF (ABS(ZLI(ISTEP)).GT.1.E-20) ZT=ZT+1./(TPCJ*ZLI(ISTEP))
  334.       IF (ABS(ZLR(ISTEP)).GT.1.E-20) ZT=ZT+1./(ZLR(ISTEP)*WLAM)
  335.       ZT=1./ZT
  336.       GO TO 16
  337. 14      ZT=CMPLX(ZLR(ISTEP),ZLI(ISTEP))/SI(I)
  338.       GO TO 16
  339. C***
  340. C***  FIX FOR PERM. OF WIRES
  341. C***
  342. C15     CONTINUE
  343. C15      RMU = ZLI(STEP)
  344. 15      RMU = ZLI(ISTEP)
  345.         IF(RMU.EQ.0)RMU = 1.
  346. C       ZT=ZINT(ZLR(ISTEP)*WLAM,BI(I))
  347.         ZT=ZINT(ZLR(ISTEP)*WLAM,RMU,BI(I))
  348. C***
  349. C***  END FIX
  350. C***
  351. 16      CONTINUE
  352.       ZREAL=ABS(DREAL(ZARRAY(I)))
  353.       ZIMAG=ABS(DIMAG(ZARRAY(I)))
  354.       IF ((ZREAL+ZIMAG).GT.1.D-20) IWARN=1
  355.       ZARRAY(I)=ZARRAY(I)+ZT
  356. 17      CONTINUE
  357.       IF (ICHK.NE.0) GO TO 18
  358.       WRITE(IW,28)  LDTAGS
  359.       STOP
  360. C
  361. C     PRINTING THE SEGMENT LOADING DATA, JUMP TO PROPER PRINT
  362. C
  363. 18      GOTO(19,20,21,22,23,24), JUMP
  364. 19      CONTINUE
  365.       ADUM=' SERIES             '
  366.       IDUM=2
  367.       CALL PRNT(ZLR(ISTEP),ZLI(ISTEP),ZLC(ISTEP),FZERO,FZERO,FZERO,
  368.      1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
  369.       GO TO 2
  370. 20      CONTINUE
  371.       ADUM='PARALLEL            '
  372.       IDUM=2
  373.       CALL PRNT(ZLR(ISTEP),ZLI(ISTEP),ZLC(ISTEP),FZERO,FZERO,FZERO,
  374.      1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
  375.       GOTO 2
  376. 21      CONTINUE
  377.       ADUM='SERIES (PER METER)  '
  378.       IDUM=5
  379.       CALL PRNT(ZLR(ISTEP),ZLI(ISTEP),ZLC(ISTEP),FZERO,FZERO,FZERO,
  380.      1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
  381.       GOTO 2
  382. 22      CONTINUE
  383.       ADUM='PARALLEL (PER METER)'
  384.       IDUM=5
  385.       CALL PRNT(ZLR(ISTEP),ZLI(ISTEP),ZLC(ISTEP),FZERO,FZERO,FZERO,
  386.      1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
  387.       GOTO 2
  388. 23      CONTINUE
  389.       ADUM='FIXED IMPEDANCE     '
  390.       IDUM=4
  391.       CALL PRNT(FZERO,FZERO,FZERO,ZLR(ISTEP),ZLI(ISTEP),FZERO,
  392.      1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
  393.       GOTO 2
  394. 24      CONTINUE
  395.       ADUM='  WIRE              '
  396.       IDUM=2
  397.       CALL PRNT(FZERO,FZERO,FZERO,FZERO,FZERO,ZLR(ISTEP),
  398.      1 LDTAGS,LDTAGF(ISTEP),LDTAGT(ISTEP),IDUM,IW,ADUM)
  399.       GOTO 2
  400. C
  401. 25      FORMAT(//,7X,'LOCATION',10X,'RESISTANCE',3X,'INDUCTANCE',2X,
  402.      1 'CAPACITANCE',7X,'IMPEDANCE (OHMS)',5X,'CONDUCTIVITY',4X,
  403.      2 'TYPE',/,4X,'ITAG',' FROM THRU',10X,'OHMS',8X,'HENRYS',7X,
  404.      3 'FARADS',8X,'REAL',6X,'IMAGINARY',4X,'MHOS/METER')
  405. 26      FORMAT(/,10X,'NOTE, SOME OF THE ABOVE SEGMENTS HAVE BEEN LOADED',
  406.      1' TWICE - IMPEDANCES ADDED')
  407. 27      FORMAT(/,10X,'IMPROPER LOAD TYPE CHOOSEN, REQUESTED TYPE IS ',
  408.      1 I3)
  409. 28      FORMAT(/,10X,'LOADING DATA CARD ERROR, NO SEGMENT HAS AN ITAG ',
  410.      1 I5)
  411. 29      FORMAT(' ERROR - LOADING MAY NOT BE ADDED TO SEGMENTS IN N.G.F.',
  412.      1 ' SECTION')
  413.       END
  414. C
  415. C
  416. C
  417.       SUBROUTINE PRNT(FL1,FL2,FL3,FL4,FL5,FL6,IN1,IN2,IN3,ICHAR,
  418.      1 IW,IA)
  419. C
  420. C     PRNT SETS UP THE PRINT FORMATS FOR IMPEDANCE LOADING
  421. C
  422.       CHARACTER ADUM*20,HALL*4,IFORM*6,IVAR*6,FMT*78,IA*4
  423.       INTEGER*4 INTT
  424.       DIMENSION IVAR(13),IA(5),IFORM(8),IN(3),INTT(3),FL(6),FLT(6)
  425.       EQUIVALENCE (FMT,IVAR(1))
  426.       DATA IFORM/'(/3X,','I5,','5X,','A5,','E13.4,','13X,','3X,',
  427.      1 '5A4)'/
  428. C
  429. C     NUMBER OF CHARACTERS PER REAL*4 VARIABLE IS 4
  430. C
  431.       DATA HALL/' ALL'/
  432. C**
  433. C     D      WRITE(*,*) '   PRNT: START'
  434. C**
  435.       IN(1)=IN1
  436.       IN(2)=IN2
  437.       IN(3)=IN3
  438.       FL(1)=FL1
  439.       FL(2)=FL2
  440.       FL(3)=FL3
  441.       FL(4)=FL4
  442.       FL(5)=FL5
  443.       FL(6)=FL6
  444. C
  445. C     INTEGER FORMAT
  446. C
  447.       NINT=0
  448.       IVAR(1)=IFORM(1)
  449.       K=1
  450.       I1=1
  451.       IF (.NOT.(IN1.EQ.0.AND.IN2.EQ.0.AND.IN3.EQ.0)) GO TO 1
  452.       INTT(1)=HALL
  453.       NINT=1
  454.       I1=2
  455.       K=K+1
  456.       IVAR(K)=IFORM(4)
  457. 1     DO 3 I=I1,3
  458.       K=K+1
  459.       IF (IN(I).EQ.0) GO TO 2
  460.       NINT=NINT+1
  461.       INTT(NINT)=IN(I)
  462.       IVAR(K)=IFORM(2)
  463.       GO TO 3
  464. 2     IVAR(K)=IFORM(3)
  465. 3     CONTINUE
  466.       K=K+1
  467.       IVAR(K)=IFORM(7)
  468. C
  469. C     FLOATING POINT FORMAT
  470. C
  471.       NFLT=0
  472.       DO 5 I=1,6
  473.       K=K+1
  474.       IF (ABS(FL(I)).LT.1.E-20) GO TO 4
  475.       NFLT=NFLT+1
  476.       FLT(NFLT)=FL(I)
  477.       IVAR(K)=IFORM(5)
  478.       GO TO 5
  479. 4     IVAR(K)=IFORM(6)
  480. 5     CONTINUE
  481.       K=K+1
  482.       IVAR(K)=IFORM(7)
  483.       K=K+1
  484.       IVAR(K)=IFORM(8)
  485.       WRITE(IW,FMT) (INTT(I),I=1,NINT),(FLT(J),J=1,NFLT),(IA(L),L=1,
  486.      1 ICHAR)
  487. C**
  488. C     D      WRITE(*,*) '   PRNT: RETURN  IW=',IW
  489. C**
  490.       RETURN
  491.       END
  492. C
  493. C
  494. C
  495.       SUBROUTINE GFOUT(CM,ZARRAY,X,Y,Z,SI,BI,ALP,BET,SALP,
  496.      1 ICON1,ICON2,ITAG,IP,IW,IGFL,LD,LD2,IRESRV)
  497. C
  498. C     WRITE N.G.F. FILE TO 'TAPE.[IGFL]'
  499. C
  500. CLARGE: CM
  501.       COMPLEX CM
  502.       COMPLEX*16 ZARRAY,SSX
  503.       COMPLEX*16 ZRATI,ZRATI2,T1,FRATI
  504. C**
  505.       REAL*4 DXA,DYA,XSA,YSA
  506.       COMPLEX*8 AR1,AR2,AR3,EPSCF
  507. C**
  508.       INTEGER*4 ICON1,ICON2,ITAG,N1,N2,N,NP,M1,M2,M,MP,IPSYM,
  509.      1 IOUT,IDM1,NEQ,NOP,NPEQ
  510.       INTEGER*4 IMAT,NPBLK,NLAST,NLSYM,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  511.       COMMON/DATAD/N1,N2,N,NP,M1,M2,M,MP,IPSYM,WLAM
  512.       COMMON /GND/ZRATI,ZRATI2,FRATI,CL,CH,SCRWL,SCRWR,NRADL,KSYMP,
  513.      1 IFAR,IPERF,T1,T2
  514.       COMMON /GGRID/AR1(11,10,4),AR2(17,5,4),AR3(9,8,4),EPSCF,DXA(3),
  515.      1 DYA(3),XSA(3),YSA(3),NXA(3),NYA(3)
  516.       COMMON/MATPAR/ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT,
  517.      1 ICASX,NBBX,NPBX,NLBX,NBBL,NPBL,NLBL
  518.       COMMON /SMAT/ SSX(16,16)
  519.       COMMON /ZLOAD/ NLOAD,NLODF
  520.       COMMON /SAVE/ KCOM,COM(19,5),EPSR,SIG,SCRWLT,SCRWRT,FMHZ
  521.       DIMENSION CM(IRESRV),ZARRAY(LD)
  522.       DIMENSION X(LD),Y(LD),Z(LD),SI(LD),BI(LD),ALP(LD),BET(LD),
  523.      1 SALP(LD),ICON1(LD),ICON2(LD),ITAG(LD),IP(LD2)
  524. C**      DATA IGFL/20/
  525. C**
  526. C $NODEBUG
  527. C**
  528. C     D      WRITE(*,*) '   GFOUT: START'
  529. C**
  530.       NEQ=N+2*M
  531. C**
  532. $DEBUG
  533. C**
  534.       NPEQ=NP+2*MP
  535.       NOP=NEQ/NPEQ
  536.       IDM1=1
  537. 50      CONTINUE
  538. C**      OPEN(IGFL,FILE='TAPE.20',FORM='UNFORMATTED',STATUS='UNKNOWN',
  539. C**  1 ERR=100)
  540. 100      CONTINUE
  541.       WRITE(*,'(A,I2,A)') ' OPEN UNIT ',IGFL,' FOR N.G.F. OUTPUT FILE'
  542. 200      CONTINUE
  543. C**
  544.       WRITE (IGFL) N,NP,M,MP,WLAM,FMHZ,IPSYM,KSYMP,IPERF,NRADL,EPSR,
  545.      1SIG,SCRWLT,SCRWRT,NLOAD,KCOM
  546.       IF (N.EQ.0) GO TO 1
  547.       WRITE (IGFL) (X(I),I=1,N),(Y(I),I=1,N),(Z(I),I=1,N)
  548.       WRITE (IGFL) (SI(I),I=1,N),(BI(I),I=1,N),(ALP(I),I=1,N)
  549.       WRITE (IGFL) (BET(I),I=1,N),(SALP(I),I=1,N)
  550.       WRITE (IGFL) (ICON1(I),I=1,N),(ICON2(I),I=1,N)
  551.       WRITE (IGFL) (ITAG(I),I=1,N)
  552.       IF (NLOAD.GT.0) WRITE (IGFL) (ZARRAY(I),I=1,N)
  553. 1     IF (M.EQ.0) GO TO 2
  554.       J=LD-M+1
  555.       WRITE (IGFL) (X(I),I=J,LD),(Y(I),I=J,LD),(Z(I),I=J,LD)
  556.       WRITE (IGFL) (SI(I),I=J,LD),(BI(I),I=J,LD),(ALP(I),I=J,LD)
  557.       WRITE (IGFL) (BET(I),I=J,LD),(SALP(I),I=J,LD)
  558.       WRITE (IGFL) (ICON1(I),I=J,LD),(ICON2(I),I=J,LD)
  559.       WRITE (IGFL) (ITAG(I),I=J,LD)
  560. 2     WRITE (IGFL) ICASE,NBLOKS,NPBLK,NLAST,NBLSYM,NPSYM,NLSYM,IMAT
  561.       IF (IPERF.EQ.2) WRITE (IGFL) AR1,AR2,AR3,EPSCF,DXA,DYA,XSA,YSA,NXA
  562.      1,NYA
  563.       IF (NOP.GT.1) WRITE (IGFL) ((SSX(I,J),I=1,NOP),J=1,NOP)
  564.       WRITE (IGFL) (IP(I),I=1,NEQ),COM
  565.       IF (ICASE.GT.2) GO TO 3
  566.       IOUT=NEQ*NPEQ
  567.       WRITE (IGFL) (CM(I),I=1,IOUT)
  568.       GO TO 12
  569. 3     IF (ICASE.NE.4) GO TO 5
  570.       REWIND 13
  571.       I=NPEQ*NPEQ
  572.       DO 4 K=1,NOP
  573.       READ (13) (CM(J),J=1,I)
  574. 4     WRITE (IGFL) (CM(J),J=1,I)
  575.       REWIND 13
  576.       GO TO 12
  577. 5     REWIND 13
  578.       REWIND 14
  579.       IF (ICASE.EQ.5) GO TO 8
  580.       IOUT=NPBLK*NEQ*2
  581.       DO 6 I=1,NBLOKS
  582.       CALL BLCKIN (CM,IDM1,IOUT,1,201,13)
  583. 6     CALL BLCKOT (CM,IDM1,IOUT,1,202,IGFL)
  584.       DO 7 I=1,NBLOKS
  585.       CALL BLCKIN (CM,IDM1,IOUT,1,203,14)
  586. 7     CALL BLCKOT (CM,IDM1,IOUT,1,204,IGFL)
  587.       GO TO 12
  588. 8     IOUT=NPSYM*NPEQ*2
  589.       DO 11 IOP=1,NOP
  590.       DO 9 I=1,NBLSYM
  591.       CALL BLCKIN (CM,IDM1,IOUT,1,205,13)
  592. 9     CALL BLCKOT (CM,IDM1,IOUT,1,206,IGFL)
  593.       DO 10 I=1,NBLSYM
  594.       CALL BLCKIN (CM,IDM1,IOUT,1,207,14)
  595. 10    CALL BLCKOT (CM,IDM1,IOUT,1,208,IGFL)
  596. 11    CONTINUE
  597.       REWIND 13
  598.       REWIND 14
  599. 12    REWIND IGFL
  600.       WRITE(IW,13)  IGFL,IMAT
  601.       WRITE(* ,13)  IGFL,IMAT
  602. C**
  603.       CLOSE(IGFL)
  604. C**
  605. C     D      WRITE(*,*) '   GFOUT: CLOSE N.G.F. OUTPUT FILE AND RETURN'
  606. C**
  607.       RETURN
  608. C
  609. 13    FORMAT(///,' **** NUMERICAL GREEN''S FUNCTION FILE ON TAPE.',I2,
  610.      1' ****',/,8X,'MATRIX STORAGE -',I7,' COMPLEX NUMBERS',///)
  611.       END
  612. C
  613. C
  614. C
  615.       SUBROUTINE COUPLE (IW,CUR,WLAM,LD,LD3,ITAG)
  616. C
  617. C     COUPLE COMPUTES THE MAXIMUM COUPLING BETWEEN PAIRS OF SEGMENTS.
  618. C
  619.       INTEGER*4 ITAG
  620.       REAL*8 DB10,GMAX,DBC,C
  621. CLARGE: CUR
  622.       COMPLEX CUR
  623.       COMPLEX*16 Y11A,Y12A,Y11,Y12,Y22,YL,YIN,ZL,ZIN,RHO
  624.       COMPLEX*16 VQD,VSANT,VQDS
  625.       COMMON /YPARM/ NCOUP,ICOUP,NCTAG(5),NCSEG(5),Y11A(5),Y12A(20)
  626.       COMMON /VSORC/ VQD(30),VSANT(30),VQDS(30),IVQD(30),ISANT(30),IQDS(
  627.      130),NVQD,NSANT,NQDS
  628.       DIMENSION CUR(LD3),ITAG(LD)
  629. C**
  630. C     D      WRITE(*,*) '  COUPLE: START'
  631. C**
  632.       IF (NSANT.NE.1.OR.NVQD.NE.0) RETURN
  633.       J=ISEGNO(NCTAG(ICOUP+1),NCSEG(ICOUP+1),LD,ITAG)
  634.       IF (J.NE.ISANT(1)) RETURN
  635.       ICOUP=ICOUP+1
  636.       ZIN=VSANT(1)
  637.       Y11A(ICOUP)=CUR(J)*WLAM/ZIN
  638.       L1=(ICOUP-1)*(NCOUP-1)
  639.       DO 1 I=1,NCOUP
  640.       IF (I.EQ.ICOUP) GO TO 1
  641.       K=ISEGNO(NCTAG(I),NCSEG(I),LD,ITAG)
  642.       L1=L1+1
  643.       Y12A(L1)=CUR(K)*WLAM/ZIN
  644. 1     CONTINUE
  645.       IF (ICOUP.LT.NCOUP) RETURN
  646.       WRITE(IW,6)
  647.       NPM1=NCOUP-1
  648.       DO 5 I=1,NPM1
  649.       ITT1=NCTAG(I)
  650.       ITS1=NCSEG(I)
  651.       ISG1=ISEGNO(ITT1,ITS1,LD,ITAG)
  652.       L1=I+1
  653.       DO 5 J=L1,NCOUP
  654.       ITT2=NCTAG(J)
  655.       ITS2=NCSEG(J)
  656.       ISG2=ISEGNO(ITT2,ITS2,LD,ITAG)
  657.       J1=J+(I-1)*NPM1-1
  658.       J2=I+(J-1)*NPM1
  659.       Y11=Y11A(I)
  660.       Y22=Y11A(J)
  661.       Y12=.5*(Y12A(J1)+Y12A(J2))
  662.       YIN=Y12*Y12
  663. C      DBC=CABS(YIN)
  664.       DBC=ZABS(YIN)
  665.       C=DBC/(2.*DREAL(Y11)*DREAL(Y22)-DREAL(YIN))
  666.       IF (C.LT.0..OR.C.GT.1.) GO TO 4
  667.       IF (C.LT..01) GO TO 2
  668.       GMAX=(1.-DSQRT(1.-C*C))/C
  669.       GO TO 3
  670. 2     GMAX=.5*(C+.25*C*C*C)
  671. 3     RHO=GMAX*DCONJG(YIN)/DBC
  672.       YL=((1.-RHO)/(1.+RHO)+1.)*DREAL(Y22)-Y22
  673.       ZL=1./YL
  674.       YIN=Y11-YIN/(Y22+YL)
  675.       ZIN=1./YIN
  676.       DBC=DB10(GMAX)
  677.       WRITE(IW,7)  ITT1,ITS1,ISG1,ITT2,ITS2,ISG2,DBC,ZL,ZIN
  678.       GO TO 5
  679. 4     WRITE(IW,8)  ITT1,ITS1,ISG1,ITT2,ITS2,ISG2,C
  680. 5     CONTINUE
  681. C**
  682. C     D      WRITE(*,*) '  COUPLE: RETURN AT END'
  683. C**
  684.       RETURN
  685. C
  686. 6     FORMAT (///,36X,'- - - ISOLATION DATA - - -',//,6X,
  687.      1'- - COUPLING BETWEEN - -',8X,'MAXIMUM',15X,'- - - FOR MAXIMUM',
  688.      2' COUPLING - - -',/,12X,'SEG.',14X,'SEG.',3X,'COUPLING',4X,
  689.      3'LOAD IMPEDANCE (2ND SEG.)',7X,'INPUT IMPEDANCE',/,2X,'TAG/SEG.',
  690.      43X,'NO.',4X,'TAG/SEG.',3X,'NO.',6X,'(DB)',8X,'REAL',9X,'IMAG.',
  691.      5 9X,'REAL',9X,'IMAG.')
  692. 7     FORMAT (2(1X,I4,1X,I4,1X,I5,2X),F9.3,2X,1P,2(2X,E12.5,1X,E12.5))
  693. 8     FORMAT (2(1X,I4,1X,I4,1X,I5,2X),'**ERROR** COUPLING IS NOT ',
  694.      1'BETWEEN 0 AND 1. (=',1P,E12.5,1H))
  695.       END
  696. C
  697. C
  698. C
  699. C***
  700. C***  FIX TO ADD PERMEABILITY TO WIRES - 12 DEC 88 - R W ADLER
  701. C***  3 LINES CHANGED
  702. C***
  703. C     FUNCTION ZINT(SIGL,ROLAM)
  704.       FUNCTION ZINT(SIGL,RMU,ROLAM)
  705. C
  706. C     ZINT COMPUTES THE INTERNAL IMPEDANCE OF A CIRCULAR WIRE
  707. C
  708. C
  709.       REAL*8 PI,POT,TP,TPCMW,CNX,BER,BEI
  710.       COMPLEX*16 CN,ZINT
  711.       COMPLEX FJ
  712.       COMPLEX*16 CC1,CC2,CC3,CC4,CC5,CC6,CC7,CC8,CC9,CC10,CC11,CC12,
  713.      1 CC13,CC14,TH,PH,F,G,BR1,BR2
  714.       DIMENSION FJX(2), CNX(2), CCN(28)
  715.       EQUIVALENCE (FJ,FJX),(CN,CNX)
  716.       EQUIVALENCE (CC1,CCN(1)),(CC2,CCN(3)),(CC3,C
  717.      1CN(5)),(CC4,CCN(7)),(CC5,CCN(9)),(CC6,CCN(11)),(CC7,CCN(13)),
  718.      2(CC8,CCN(15)),(CC9,CCN(17)),(CC10,CCN(19)),(CC11,CCN(21)),(CC1
  719.      32,CCN(23)),(CC13,CCN(25)),(CC14,CCN(27))
  720.       DATA PI,POT,TP,TPCMU/3.1415926D0,1.5707963D0,6.2831853D0,
  721.      1 2.368705D+3/
  722.       DATA CMOTP/60.00/,FJX/0.,1./,CNX/.70710678D0,.70710678D0/
  723.       DATA CCN/6.E-7,1.9E-6,-3.4E-6,5.1E-6,-2.52E-5,0.,-9.06E-5,-9.01E-5
  724.      1,0.,-9.765E-4,.0110486,-.0110485,0.,-.3926991,1.6E-6,-3.2E-6,1.17E
  725.      2-5,-2.4E-6,3.46E-5,3.38E-5,5.E-7,2.452E-4,-1.3813E-3,1.3811E-3,-6.
  726.      325001E-2,-1.E-7,.7071068,.7071068/
  727.       TH(D)=(((((CC1*D+CC2)*D+CC3)*D+CC4)*D+CC5)*D+CC6)*D+CC7
  728.       PH(D)=(((((CC8*D+CC9)*D+CC10)*D+CC11)*D+CC12)*D+CC13)*D+CC14
  729. C      F(D)=DSQRT(POT/D)*CEXP(-CN*D+TH(-8./X))
  730.       F(D)=DSQRT(POT/D)*CDEXP(-CN*D+TH(-8./X))
  731.       G(D)=CDEXP(CN*D+TH(8./X))/DSQRT(TP*D)
  732. C***
  733. C***  FIX FOR REL. PERM.
  734. C***
  735. C     X=DSQRT(TPCMU*SIGL)*ROLAM
  736. C      X=DSQRT(TPCMU*RMU*SIGL)*ROLAM
  737.       X=SQRT(TPCMU*RMU*SIGL)*ROLAM
  738. C***
  739. C***  END FIX
  740. C***
  741.       IF (X.GT.110.) GO TO 2
  742.       IF (X.GT.8.) GO TO 1
  743.       Y=X/8.
  744.       Y=Y*Y
  745.       S=Y*Y
  746.       BER=((((((-9.01E-6*S+1.22552E-3)*S-.08349609)*S+2.6419140)*S-32.36
  747.      13456)*S+113.77778)*S-64.)*S+1.
  748.       BEI=((((((1.1346E-4*S-.01103667)*S+.52185615)*S-10.567658)*S+72.81
  749.      17777)*S-113.77778)*S+16.)*Y
  750.       BR1=DCMPLX(BER,BEI)
  751.       BER=(((((((-3.94E-6*S+4.5957E-4)*S-.02609253)*S+.66047849)*S-6.068
  752.      11481)*S+14.222222)*S-4.)*Y)*X
  753.       BEI=((((((4.609E-5*S-3.79386E-3)*S+.14677204)*S-2.3116751)*S+11.37
  754.      17778)*S-10.666667)*S+.5)*X
  755.       BR2=DCMPLX(BER,BEI)
  756.       BR1=BR1/BR2
  757.       GO TO 3
  758. 1     BR2=FJ*F(X)/PI
  759.       BR1=G(X)+BR2
  760.       BR2=G(X)*PH(8./X)-BR2*PH(-8./X)
  761.       BR1=BR1/BR2
  762.       GO TO 3
  763. 2     BR1=DCMPLX(.70710678,-.70710678)
  764. C***
  765. C***  FIX FOR REL. PERM.
  766. C***
  767. C3     ZINT=FJ*DSQRT(CMOTP/SIGL)*BR1/ROLAM
  768. C3     ZINT=FJ*DSQRT(CMOTP*RMU/SIGL)*BR1/ROLAM
  769. 3     ZINT=FJ*SQRT(CMOTP*RMU/SIGL)*BR1/ROLAM
  770. C***
  771. C***  END FIX
  772. C***
  773.       RETURN
  774.       END
  775.